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We propose a simple model for a binary decision making process on a graph, motivated 
by modeling social decision making with cooperative individuals. The model is similar to 
a random field Ising model or fiber bundle model, but with key differences in behavior on 
heterogeneous networks. For many types of disorder and interactions between the nodes, 
we predict with mean field theory discontinuous phase transitions which are largely inde- 
pendent of network structure. We show how these phase transitions can also be under- 
stood by studying microscopic avalanches, and describe how network structure enhances 
fluctuations in the distribution of avalanches. We suggest theoretically the existence of a 
"glassy" spectrum of equilibria associated with a typical phase, even on infinite graphs, 
so long as the first moment of the degree distribution is finite. This behavior implies that 
the model is robust against noise below a certain scale, and also that phase transitions 
can switch from discontinuous to continuous on networks with too few edges. Numerical 
simulations suggest that our theory is accurate. 
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1. Introduction 



Over the past decade there has been an explosion of interest in the statistical physics community into the 
behavior of simple statistical models on networks. These models often lead to insights about qualitative 
behavior in social systems, as random graphs are a low order approximation to realistic social networks 
[31 [7] . Decision making processes have long been studied as a simple example of such an application. 
The voter model [27] or variations with nonlinearities [9] or other complications [30, 10J have long been 
proposed, although they are only models of consensus building. The Axelrod model is an alternative 
model which instead exhibits equilibria with diversity of opinions J2[ El [2"5j |14j . Other spin models or 
agent-based models have been proposed to study financial markets [51 [25] . 

Many of the above models do not predict a key phenomenon, however, which is the presence of 
shocks, catastrophes and discontinuous phase transitions as external parameters are slowly tuned. Often 
times, entirely new models have been proposed simply to account for this phenomenon |31| |23] 111] . 
More interesting, however, is the proposal that a random field Ising model, well-known for hysteresis and 
discontinuous phase transitions [26], can be used to model these phenomena in social science [U [61 [T3] . 
A similar model called the fiber bundle model, used to study the breakdown of some materials, also has 
similarly promising features [24, 16J. It has been noted [29J that there are some similarities between the 
fiber bundle model at the critical point and financial market failures, although beyond the critical point 
the analogy has not been studied significantly. Other recent models have attempted to discuss disorder- 
induced phase transitions of opinion dynamics, using disorder in the interactions between individuals 



In this paper, we propose from scratch a very simple model for a binary decision making process on a 
network, in which individuals make reactionary decisions as a global external field is changed. Our result 
will somewhat look like a generalized fiber bundle model, or alternatively a random field Ising model 
in a global magnetic field, but with some important differences which make our model nearly exactly 
solvable on heterogeneous graphs. We provide a preliminary mean field analysis of the model predicts 
disorder-induced discontinuous phase transitions and hysteresis, and provide a microscopic justification 
for mean field theory as well, describing the microscopic dynamics of binary decision making in terms of 
avalanches. We show how heterogeneous networks with heavy-tailed distributions can cause the avalanche 
dynamics to have fluctuations which scale with the size of the network. 

Most interestingly, however, we will show that there is an infinite spectrum of equilibria (in the large 
graph limit). Again, this is not too surprising, because the random field Ising model has spin glass- 
like characteristics on networks which are more enhanced on scale free networks |17| [2"T] [T5] (although 
mathematically by a standard definition, it is not a spin glass [20J). Since our model does not admit a 
Hamiltonian and therefore a free energy, the "glassy" behavior of the binary decision model will instead be 
characterized by the presence of this spectrum of equilibria. We will then use this spectrum of equilibria to 
justify two very interesting phenomena: the robustness of the binary decision model to small fluctuations, 
and the possibility that network structure can suppress a discontinuity in the phase transition. The 
realization and quantitative predictions of "glassy" multistability proposed in our binary decision model 
is new in the literature, to the best of our knowledge. 
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2. The Binary Decision Model 



In this section, we will introduce the binary decision model, justifying its use as a simple model for 
equilibrium social behavior. We begin by approximating that individuals interact via a social network, 
which can be described as an undirected graph G consisting of a vertex set V and edges E. We will 
denote N = \ V\: i.e., there are N nodes in the graph. To each node v in the graph, we associate a binary 
variable x v G {0, 1}. For example, x v = may mean that the individual is uninterested in participating 
and trading in a given economic sector, while x v = 1 means the opposite; alternatively x v = could 
model that an individual does not have an active account for a social media service, with x v = 1 the 
opposite. Each node will decide its state, or 1, by comparing an internal field, which we label s v , to an 
external (global) field p^ 

x v = @(s v -p). (1) 

For example, if we think of p as the external price of some good, then s v represents the effective price at 
which buyer v is willing to buy: when s v > p, x v = 1 and the buyer is actively buying, and when s v < p, 
x v = and the buyer is not actively buying. 

In order to fully specify the model we thus simply need to describe how to determine s v . Our 
formulation of the binary decision model will be to assume that 

s v = P v h(q v ) (2) 

where P v is some internal variable, h is a monotonically increasing function, and 

q v = P(x u = l\uv € E) (3) 

simply represents the fraction of neighbors of v which are in state 1. We are free to rescale P v so that 
h(0) = 1, and we will assume so for the remainder of the paper. Note that we will always assume that h, 
P v and p are non- negative. Note that under the identifications 

p = eP+, (4a) 
Pv = e P+ ", (4b) 
h(q)=e h+iq \ (4c) 

we can formally write the binary decision model as 

x v = @(h + (q v ) + P +v -p+), (5) 

since the exponential is a monotonically increasing function. 

To understand the choice above, it is helpful to re-cast the problem temporarily in the language 
of economics. Suppose that p represents some globally observed price (or, more abstractly, "utility" 
parameter); then, P v represents the price/utility that node v believes is the value of existing in state 1. 
The key point of this model is that the effective value of P v is multiplied when neighbors of v are also 
in state 1. For example, a social media service (where P represents a utility, not a specific price) is far 
more valuable to its users when there are many other users. Similarly, stock traders may base much of 
their valuation of a stock based on what they believe the rest of the market is doing. Therefore, this 
binary decision model represents social scenarios where individuals are making a binary decision based on 
comparing the utility of two options, when that utility is dependent on how their neighbors are behaving. 



1 It is not important what happens when s v — p for almost all reasonable formulations of the graph G, or the internal 
fields s v . 
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The argument can and has been made [6] that an alternative choice, 

s v = P v + J x u , (6) 

(uv)eE 

which directly results in the famous random field Ising model, is also worth studying. In fact, note that 
this model is formally equivalent to the binary decision model on graphs where all nodes have the same 
number of edges, for a special choice of h+(q) = kJq. For the purposes of this paper, we will stick with 
([2]), which leads to simpler calculations. We should note, however, that there is one major drawback to 
a choice such as this: as far as we can tell, there is, in general, no Hamiltonian function which will result 
in the binary decision model, and this denies us the use of some of the tools of statistical mechanics. 



3. Mean Field Theory 



Let us now turn to mean field theory to "solve" the binary decision model, and compare to simulations. 
Our first pass will demonstrate both how mean field theory can be quite accurate, as well as reveal some 
of its major drawbacks. 



3.1. A Macroscopic Solution 

Let us denote 



P(x„ = 1) (7) 



where the average over nodes is over the uniform distribution over V. Let us also assume that the random 
variables P v are i.i.d. drawn from a probability distribution with cumulative distribution function 1— P(P), 
and probability density function /(P). It is straightforward to see that 

q = P(h(q)P v >p)=F (JL^J (8) 

where the mean field approximation that q v = q has been applied. Q is the mean field equation describing 
the possible equilibria of the system, as described by the single parameter q. Note that this mean field 
equation appears independent of the degree distribution of the graph, if we chose to take this into account, 
since 

q k = P(x v = l\k v = k)= P(x v = l)=q. (9) 

The middle equality above is a consequence of the fact that in mean field theory, every edge points to 
the same "mean field node": i.e., every edge contributes the same factor of q. There are cases where this 
approximation will break down, and we will return to this at the end of the paper. 

We are most interested in the case where there are multiple equilibria, i.e. where Q has multiple 
solutions. Graphically it is clear how to approach this problem, as shown in Figure [T] Inspired by this 
approach, we can also straightforwardly analyze it analytically. Since F is a CDF, we have < F < 1. 
Because we must have < F(q = 0) and 1 > F{q = 1), there must be some q* for which Q is true. 
Furthermore, let us consider the behavior of F near this point. Suppose that dF/dq > 1 - this means 
that F < q for q just smaller than q* , and F > q for q just larger than q* . However, given the constraints 
at q = and q = 1, we see that there must be at least 2 other points at which q = F, implying the 
existence of at least 3 solutions to (|8l) . We conclude that if 



dF(p/h(q)) 



dq 



„. =f yw = ° i <io) 
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one solution: critical point: bistability: 




0.5 1 0.5 1 0.5 1 



q q q 

Figure 1: The graphical method of solving equation Q. The solid dots represent equilibria which 
are stable , and the empty dots represent equilibria which are unstable: the half-filled dot represents a 
marginal point. 



then we have multiple equilibria. In a later section, we will find a simple microscopic interpretation for a 
as the probability that a spontaneous flip in any node's state will cause one of its neighbors to also flip, 
and the stability criterion a < 1 follows from the condition that an avalanche have finite expected size. 

Now, let us consider what happens as we change p. In particular, suppose we increase p to some 
p = p c for which 

h(q*)J h(q*) 2 1 ' 

Then we conclude that passing through p c in the direction which decreases the left hand side of ( |11[ ) 
corresponds to the disappearance of a pair of fixed points. In particular, if the solution at q* disappears, 
it must discontinuously jump to another point. This is the hallmark of a discontinuous phase transition. It 
is not always the case that a distribution of F(P), and an interaction h(q), will allow such a discontinuous 
phase transition, but quite often discontinuous phase transitions do occur. Essentially, increasing the 
value of h(q) makes the model more and more likely to admit a discontinuous phase transition, and 
subsequently increase the size of the jump at the transition. 



3.2. An Exactly Solvable Case 

Let us look at a sample of an exactly solvable version of this model, to the level of approximation we have 
just studied. We will take 

h(q) = l + Aq. (12) 

and 

r i p<p 

Po + l-P P <P<P + 1 . (13) 
^ P>P + l 

It is easy to check when q = is a solution: this occurs when F(p/h(0)) = 0, or F(p) = 0, or p > Pq + 1. 
q = 1 is a solution if F(p/ (1 + A)) = 1, or p < (1 + A)Pq. For < q < 1, solutions occur when 



F(P) 



1 + Po 



P 



l + Aq 

which can easily be solved using the quadratic formula to give 



^o + l 



1 



Po + l + 



A 



4p 
~A 



(14) 



(15) 
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The important feature of (15) is the square root, which will become imaginary at price 

2 



Pc 



A 

7 



(16) 



There are 3 distinct qualitative possibilities for the phase diagram. To understand these possibilities, 
it will suffice to consider q at p = p c , which is 



Po + l 



(17) 



If g c > 1, which occurs when Pq > 1 + then we know that only one branch of physical solutions exists 
(other than q = or 1), and this branch has a positive slope in the pq plane, connecting p = (1 + A)Pq 
and p = Pq + 1. If q c < 0, then one branch of physical solutions exists, with negative slope, connecting the 
same 2 points. Otherwise, we see that q c is a physical value, and therefore is a valid point in parameter 
space - in particular, there are two branches of allowed solutions (other than q = or 1). We show 
examples in Figure [2j 



P Q = 0,A = 0.5 



P = 0, A = 4 



3, A = 1 





1 - 
0.5 - 

- 



4 6 
P 



Figure 2: Three examples of exact mean field solutions for uniform distributions on P v displaying 
qualitatively different behaviors. The solid lines represent stable equilibria, and the dashed lines unstable 
equilbria - we will understand stability from a microscopic standpoint later. 

Physically, we see that depending on the parameters, the binary decision model has a variety of 
interesting behaviors. When A is small enough, there are no discontinuous phase transitions, suggesting 
"nice behavior". However, when A gets large enough, discontinuous transitions begin to occur, and when 
A reaches a critical value, the only stable equilibria are with the entire graph in state or 1, representing 
a hyper-polarized situation. 



3.3. Numerical Simulations 

We now present numerical simulations of the binary decision model. In our numerical simulations, we 



always take h(q) to have the linear form (12), although we allow the coefficient A to vary. We ran 
simulations over both Erdos-Renyi random graphs, drawn from an ensemble where every possible edge 
is equally likely to be included in the graph, as well as scale free graphs, generated by the algorithm of 
|19| . For our internal P distributions, we generated P v from either uniform distributions over [0,1], or 
Gaussian distributions with \x = 1 and a ~ 0.2 (with varying a), and scale free distributions with varying 
exponents and minimum P v of 1. To observe bistability if it exists, we began our simulation by starting 
with p = 0, and then increased p in uniform steps up until some value p maX ) and then decreased p in equal 
steps back to p = 0. 
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As a first simulation, we demonstrate that for a given internal disorder distribution F(P), and a given 
h(q), the mean field theory solution ^ is typically a very good approximation. This is shown in Figure 
[3] Although the mean field theory is typically quite good, we already notice some key deviations from 
our prediction. Note that typically, the phase transition does not look as sharp as it occurred because 
we are averaging over runs, and there are small fluctuations in the critical value of p due to internal 
disorder. Figure [4] shows the emergence of a phase transition as the interaction strength A is increased 
past a critical value of 1, for the given uniform price distribution, as shown by both theory and numerics, 
as well as some sample distributions where P v is a Gaussian or scale free random variable. 



1 



Q 0.5 





0.5 1 1.5 2 

P 

Figure 3: The average value of q(p) for a uniform distribution of prices with A = 2. We averaged over at 
least 50 trials for each type of network, and ran over both Erdos-Renyi and scale free graphs with varying 
values of N, (k), and scale free exponent - each graph ensemble is represented by a different color. The 
smoother transitions correspond to graphs with smaller (k), a fact which we will explain later. 

While the mean field theory approximation appears to typically hold, we should point out some key 
failures. In particular, there is a noticeable hysteresis effect associated to the upper branch of the mean 
field solution near the phase transition. We also notice that the phase transition appears to be slightly 
delayed for larger values of k, and on some graphs it appears as though there is no discontinuous phase 
transition at all. Both of these phenomena are key signatures of the advertised multistability, and we will 
return to them later. 

3.4. Finite Size Fluctuations from Disorder 

Let us briefly discuss the fluctuations about mean field theory due to finite size effects (but not small 
network effects). These fluctuations are simply due to the fluctuations in the values of the internal fields, 
P v . To quantitatively estimate their size, let us denote q = qo + 5, with qo the mean field value predicted 
by ([8]) and 5 a small fluctuation. Similarly, let us denote the CDF of the actual distribution of P, realized 
on the graph, as F = Fq + A, with Fq the mean field value and A a small fluctuation. Then, expanding 
^ to lowest order in the fluctuations we find 




which implies that 

A 
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Figure 4: Here we show how increasing the "interaction strength" in the binary decision model causes 
the onset of a discontinuous phase transition, both in theory and simulation. Our simulations for this 
graph used Erdos-Renyi graphs with 5000 nodes and {k) = 20: the distributions in the graphs refer to 
the distributions of the internal disorder. Note that the discontinuous phase transition appears smeared 
out because of fluctuations in the realizations of disorder. 



The distribution of A is simple to find, although we will focus only on the variance of the fluctuations, 
as the higher order fluctuations are rapidly suppressed for increasing N. Since the internal disorder 
consists of iid random variables, and for the given value of p and the distribution of P v the probability 
that a node is in state x v = 1 is Fq, assuming the rest of the network to be in a mean field state, we 
conclude 



Var 

This implies that 



^>)-^£v^)-^-*L=4. (20, 



Since this number is typically quite small due to the factor of 1/N, and our theory suggests that this 
variance is not dominated by large deviations, we postulate that higher order fluctuations will not play 
an important role. We find that this relation is obeyed very well so long as we do not approach a = 1, as 
shown in Figure [5| 
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Figure 5: We compare the predictions of (21) to numerical simulations on a variety of graphs. Different 
shapes represent different values of A, or different distributions F(P), some of which admit discontinuous 
phase transitions, and others which do not. Different colors represent different values of iV (ranging 
from 800 to 5000) and/or degree distribution (Erdos-Renyi and scale free with exponents 2.5 < 7 < 4). 
Significant deviations come near the onset of a phase transition - we have removed some of these values 
near the critical point when our code is averaging over realizations in different phases. 



4. Avalanches 



We now turn to a slightly different question, which is a reinterpretation of our previous mean field 
analysis on the existence of a phase transition. Let us consider a state of the graph in equilibrium at 
q* = F(po/h(q*)). Suppose that the external field po is increased slightly to p, causing only node vq to 
flip. How many other spins will also flip, due to the change of state of vq! 

Let us define Q to be the probability that one of the vertices v connected to vq will flip: 



q = p( ,, * >p«> * Wf p 



h(q v -l/kv) h(q v )J \h(q)J \ \h(q - 1/k) j , edges 

edges 



F| .^r F ^jni / ^Jw;*.- a w. (22) 



Now, we have to be a bit careful. As noted above, the averaging occurs over the distribution of nodes 
which an edge points to, which is different from the distribution of nodes itself, because nodes with more 
edges are counted more often: in fact, each node v is counted k v times in the average (• • • ) e dges over nodes 
pointed to by an edge. Therefore: 

Ocdacs = ^ljil = W (23) 
so we conclude that 

«=w- (24) 

Let us now denote with n the number of the vertices v we expect to transition to the state. For any 
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given site, this occurs with probability Q, so 



■'<) 



(n) = (J£ Qj = E ^ = <*>Q = a - ( 25 ) 

Assuming that Q stays constant, if any new spin flips, then it too will have the possibility of flipping spins. 
If we approximate the new spin as the same as the old spirj^] then, since the avalanche approximately 
grows as a birth/death process, we conclude the total size of the avalanche is 

("total) ~ -i ~ • (26) 
1 — a 

Of course, this formula only holds for a < 1; for a > 1, it is well known that the birth/death process has 
an expected infinite size, and for a = 1 it is almost surely finite but with infinite expected value. 

This gives us much new insight into the physical processes at work behind mean field theory. We see 



looking back at (11) that we found that reaching a point where a = 1 corresponds to a phase transition 
- here, that means that we expect a spin avalanche to have infinite size. Furthermore, it allows us to 
analyze the stability of the fixed points we found earlier - only the fixed points where a < 1 are stable. 
Understanding phase transitions by considering microscopic avalanche sizes is not new: see, e.g., [6] in 
the context of the random field Ising model, or |24] [16], [32], ITS] in the context of the fiber bundle model. 
Interestingly, however, for our binary decision model, this calculation helps to give insight into why the 
graph structure is seemingly so unimportant in mean field theory - even though more connected nodes 
are affected more often by spin flips, they are not affected as much, and these effects cancel each other. 
This would not be the case in the random field Ising model on a heterogeneous network, for example, 
where a itself would obtain an intricate dependence on the degree distribution. In this sense, our binary 
decision model, where h(q) is fc-independent, is a very convenient toy model to solve. 

Let us also ask about the fluctuations in the size of avalanches, which we will explore by considering 
variations in the size of n, the number of neighbors that one site flipping will also flip. These fluctuations 
occur both within the same realization of graph structure and internal disorder. However, for simplicity 
we have included many realizations in our average. Denoting z v = 1 if v flips and otherwise: 




/ k 



Var(n) = (n 2 ) - (n) 2 = ^Pk ( [ J>i ) " « 2 = Z> ( Z>J + 

v i=i / / \i=i m I 



a 2 



= Vft [kQ + k(k - l)Q 2 ] -a 2 = a(l - a) + ^—^a 2 . (27) 

We see that the graph structure thus plays an important role in fluctuations of the size of avalanches. 
This theory is quite accurate as shown in Figure |6j although it does seem to break down very close to a 
phase transition, as the diverging curves suggest. 

Other work [32] [TS] considers in more detail the theory behind avalanche distributions in the fiber bun- 
dle model, which is quite similar to the binary decision model, on fully connected networks, in particular 
near the critical point, although there is no general theory of the distribution of avalanches accounting for 
heterogeneous network structure. Indeed, our results above, as we saw that the variance included network 
structure even assuming a "mean field" network structure. 



2 Arguably, the k above should be replaced with a k — 1, but for k 2> 1 this is not a major qualitative or quantitative 
change. We also require k 2> 1 for our Taylor approximation above to be accurate, so we will for simplicity neglect worrying 
about this for this paper. 
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Figure 6: We compare theory vs. numerics for the predicted variance in the number of nodes which 
change state based on the change of a single node. Different colors represent different graph ensembles, and 
different shapes are different F(P) and/or h(q). The deviations appear to become much more significant 
as we approach a phase transition, perhaps due to loopy effects in the graph (where our avalanche theory 
breaks down). 



5. Multistability 



In this section, we will discuss the binary decision model's most interesting feature - "multistability." 
By the term multistability, we mean that there is a continuous spectrum of equilibrium states, even in 
the N — t- oo limit, so long as (k) is finite. Recall that an equilibrium for the binary decision model is 
a state where all nodes satisfy the constraint equation that x v = Q(P v h(q v ) — p). Since q v depends on 
the values of x u for each neighbor u of node v, it is possible that there are multiple possible solutions to 
the constraint equation - thus, multiple equilibria exist. We've trivially seen that this can be true by the 
existence of two phases, but here we are interested in the possibility that multiple equilibria may exist 
for any given phase. 

Due to the internal disorder, and the similarity to the random field Ising model, it is not surprising that 
glass-like behavior arises, although the lack of a formal Hamiltonian or free energy makes it challenging 
to classify the binary decision model as a glass. Instead, we will call the behavior multistability. The goal 
of this section is to propose multistability from a theoretical standpoint. We will then show it exists in 
our numerical simulations, and comment on the implications of this phenomenon. 



5.1. A Spectrum of Equilibria 

Let us approximate the probability that a pair of spins will satisfy the x v constraint either if they are 
both in state 1 or in state 0. Note that it will never be the case that we would need to consider the 
possibility that a pair of nodes could flip between 10 and 01, because the constraint equation implies 
that decisions are always made cooperatively]^] We will assume that the remainder of the graph is treated 
within the mean field approximation. The probability that one node flipping will flip another node has 



been calculated in the previous section in (24). Thus, we simply calculate the probability that two nodes 
can be in either state 1 or to be the probability node u flipping flips node v and vice versa, which is 
simply 

P(pair can exist in 2 states) = ( — - ) . (28) 

\( k )J 



3 With more work, this can be made rigorous. 
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Since there are N{k)/2 edges in the graph, we thus approximate the number of equilibria as being 



equilibria 



yV« 2 /2(A;} 



(29) 



However, we now need to take into account the expected value of the avalanche caused by the pair flip. 
This can be accounted for by roughly multiplying by a factor of (1 — a) -1 , corresponding to the expected 
size of the avalanche caused by one of the nodes. Thus, we find that given a pair of nodes, we should 
expect that flipping them will cause a change of 



Aq, 



1 



2a 2 



one pair 



N (l-a)(fc) 2 ' 



and this leads to a width of the equilibrium spectrum which is finite even in the iV 
(k) is finite: 



(30) 



oo limit, so long as 



Aq 



a 



(31) 



(l-a)(fc) 

As discussed earlier, we ran simulations of the binary decision model by first increasing p from a state 
of all 1, and then decreased p from states with many 0s. This means that we can observe, quantitatively, 
the range of the spectrum of equilibria by observing the difference between the value of q on the upward 
sweep versus the downward sweep. Figure [7] shows that (31 ) is quantitatively correct, despite the incredibly 
simple theory we used. 




10 
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Figure 7: A comparison of the numerically determined width of the equilibria spectrum to (31 ). Different 
colors represent different N (ranging from 1000 to 5000) or pk (both Erdos-Renyi and scale free with 
2.5 < 7 < 4), and different shapes represent different F(p) and h(q). We see that the theory is very 
accurate even near a phase transition (these are the points with Aq larger). 

Indeed, it should not be so surprising that we found characteristics of glassy behavior, considering the 
similarity of this model to the random field Ising model. Interestingly, it is unlikely this model is a true 
spin glass, due to its similarity to the random field Ising model where it was shown that the spin glass 
susceptibility is not divergent |20] . This is perhaps more of a mathematical technicality than a physically 
meaningful statement: in our model, so long as (k) is finite and a > 0, the binary decision model is always 
a multistable "glass" . 
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5.2. Multistability by the Cavity Method 

Since the multistability phenomenon is so fundamental to our model, let us predict multistability via a 
second approach: roughly speaking, we will derive TAP-like equations for x v via the cavity method |12j . 
The cavity method is an extension of mean field theory which is used to understand disordered systems, 
and it essentially refers to picking a special node in the graph, called the "cavity" node, and evaluating 
the probability that it is in each state, assuming that all of its neighbors feel the effect of the cavity's state 
exactly, with the remainder of the graph treated at mean field level. Of course, this can be extended: an 
m th order cavity method could take into account the reaction of the cavity on nodes up to m edges away. 
Note that the approach is also in the spirit of belief propagation on a tree [22] . 

To use the cavity method for our problem is straightforward: we need to make our mean field argument 
from earlier a bit more refined, so we will pick a special node v (the "cavity") and explicitly write 

*« = e ( p "-/^)) < 32 > 

and explicitly determine q v . To find q v , we need to determine the expected value of a state x u , for each 
u with uv G E. If we only use a first order cavity approximation, then k u — 1 of the edges of u point to 
nodes with x = q, so we find 

x u = 6 (p u - - ■ -. P -—] . (33) 



h(q + (x v - q)/ku) 

Averaging over the internal disorder and graph structure it is straightforward to show that, assuming 
that k is large, 

P(x„ = 1) = „ = F ( JL) + Jl ( „ _,)_,+ JL(„ _ ,). (34) 

Then we find the equation for the cavity node v, assuming again that k is large enough that we can Taylor 
expand h(q): 

*-K F *-^ + iwwM- (35) 

Now averaging over the internal disorder of P v , we find that it may be possible that both x v = and 
x v = 1 are solutions. The calculation of how likely this is is very simple: 

P f 1| %) tt ^ )f) P L | h'(q)a(q- _ / p /^h'(q)aq\\ a ph'(q) 



2 



h{q) V %) (k)J h(q) V h(q) (k) J J J \h(q) \ h(q) (k) J J (k) h(q) 

= ^ + 0(a 3 ) (36) 

where the 0(a 3 ) terms correspond to terms proportional to the derivative of /. This is precisely what we 
found earlier, neglecting the back reaction onto the remainder of the graph (via avalanches), up to the 
new terms which have arisen via the cavity method. 

Let us now perform the cavity method to higher orders: in particular, let us assume that the graph 
is tree-like (at least locally), and keep track of all the nodes up to a graph distance of n away from the 
cavity. The tree-like approximation is very convenient, as it allows us to assume that each node feels the 
effects of the cavity node v through exactly one neighbor. We will start by considering the case of n = 2 
- the generalization to larger n will be very straightforward. The state of the cavity is given by 



x v = 



( \ 

p.- " 



1 



h I — > x 



(37) 
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and xi, . . . , Xk v denote the states of v's neighbors. The state of one of the neighbors of v, e.g. x±, is given 
by 



xi = e 



\ 



Pi- 



p 



V 



h 



+ i 



1 



(38) 



7 



where q\ corresponds to the fraction of nodes (other than v) which are neighbors of node 1, and are in 
state 1. But note that if x\ is unknown, we already computed the formula for q\ previously: 



a 



(39) 



where as before, q is the fraction of nodes far from the graph. Now, we write, assuming the number of 
nodes is large, as usual: 



xi = e P] 



p 



h(q) 



1 



h(q) 



a 



i 



i 



(40) 



and using the same argument as before by finding the largest and smallest possible corrections to the 
effective value of p, we conclude that the probability that there is an equilibrium state with x\ = 1, as 
well as one with x\ = 0, given x v , is given by 



P(xi = or l\x v ) 



a 



a 

W) 



1 



1 



+ 0(a 2 



(41) 



Note that we have not yet allowed for fluctuations in xq, which we have assumed is fixed. Also, given 
this formula, it should be fairly clear that we could have guessed this answer a priori from what we found 
before, under the assumption that one node x v is fixed. 

Finally, let us return to the question of interest: the probability that x v can be in both states. We 
need to compute the largest and smallest possible values of q v : 



P(m 
P(m 



1 \%v 

0\x v 



l) 



q + a 



0) = q + a 



a 




Jk) 




a 




Jk) 



l-q 



Q_ 



(42a) 
(42b) 



Note that the formula above neglects some of the smaller corrections due to derivatives in /, e.g., although 
these could easily be carried through if needed. We may finally average over the degree distribution, and 
replace l/k± with l/(k). Using these equations to find upper and lower bounds on q v , we finally can 
conclude, using the same logic as in our earlier computation, that 



P(x v = or 1) = a 



a a 



1 

Jk) 



a 2 (l + a) 

W) 



(43) 



Note that this is also the width, Aq, of equilibria, up to second order in the cavity method. 

Now, let us extend this computation to higher orders. If we have found the width of the spectrum 
accounting for all nodes up to a distance n away, assuming the graph is a tree, we can easily determine 
the width of the spectrum accounting for all nodes a distance n+1 away from x v , by simply treating the 
neighbors of x v as the cavities and using the result for n. An analogous computation to the one we did 
above shows that: 



a ■ • » 



l(k) 



(44) 
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Using the results for AgW and Aq^ from before, it is clear that 



A<Z (n) = ~$> 1+fc - ( 45 ) 
W k=i 

Summing this series to infinite order, we find that 

A ? (oo) = n a \i^v ^ 

(1 -a){k) 

assuming that a < 1 so that the series is summable. Thus, we see that the cavity method recovers the 
result we found using simple logic earlier. 

5.3. Robustness to Fluctuations 

An important consequence of the spectrum of equilibria is that the macroscopic model is robust against 
small perturbations in the external field p. To estimate the size 5 of a perturbation in the external field 
p required to change the macroscopic state (the value of q)j^]we use the fact that 

dg(p) 1 " 2 U7) 
dp ~ 6(l-a)(k)' [ ' 



(47) follows directly from basic graphical considerations, by considering ([8j) and using that if the slope of 
the mean field curve is known, then the fluctuations in both the q and p directions in F(p/h(q)) must be 
related. Taking an implicit derivative of Q, we find that 

dq(p) = a h(q) 

dp 1 — a ph'(q) ' 

which implies that 

x_ "l^fo) (49) 



(k) h(q) 

For a generic model, the /i-dependent factor in 5 is likely 0(1), and so the dominant feature is the a 
and (A;) dependence. We see that the model becomes more robust against small price fluctuations as we 
approach a critical point (where a = 1), although this is a linear approximation, as we assumed that the 
fluctuations were small enough that we could neglect terms beyond the first derivative in F(p/h(q)), in 



(47). Thus, (47) will break down as nonlinear effects become important as we approach the critical point. 



5.4. Suppression of Discontinuous Phase Transitions 

What happens to the fluctuations of multistability as we approach the critical point, which mean field 
theory will in general take as a discontinuous phase transition? Here, the simple linear arguments we 
used above begin to fail. The l/(k) fluctuations can become so large on graphs with "small" values of (k) 
that they can in fact remove the discontinuity in the phase transition. 

We can use a simple mean field theoretic argument to justify this. Let q k denote the fraction of nodes 
with k edges in state 1, and r the probability that an edge points to a node in state 1: 

r = Tjry2kp k q k . (50) 



4 It is possible that the equilibrium that we are at is at an outer edge of the spectrum, where this argument does not work, 
but for most equilibria (which are in the center of the spectrum) it works fine. 
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A more realistic expression for qk is given by 



m=0 



fc! 



m\(k — m) 



-r m (l-r) m F 




(51) 



For simplicity, let us suppose that k is large enough that binomial distribution is well approximated by a 
Gaussian distribution, so we have 



Qk 



1 

/ 



dq 



v / 2^AF T r(r 



-k(q-r) 2 /2r(l-r) F 



P 

h(q) 



(52) 



Certainly when k gets very large, the Gaussian collapses to a 5 function and we obtain q^ = F(p/h(r)), 
independently of k. However, suppose that the length scale qp of fluctuations in F(p/h(q)) is comparable 
to ^fk^ril — r). Then we see above that fluctuations in the number of nodes which are actually in state 
1 may smooth out fluctuations in F, destroying the phase transition! 



Since (51) is far too complicated to analyze in the regime of interest, where k becomes small, we 



resort to estimating numerically whether or not there will be a phase transition. To do this we begin by 



determining self-consistent values of r, as calculated by using (50) and (51). It is simplest to imagine 
doing this by plotting the function r, and then the function determined by (50) and (51), and looking for 
intersections. The reason this is the preferred method is because we have not yet taken into account the 
l/(k) fluctuations. A crude way of doing this is simply adding l/2(k) to the right hand side of ( |50| ), as we 
can consider the multistability as being caused by fluctuations in the value of F(p/h(q)) of ±a/2{k). As 
the previous sections have suggested, this effect can become quite serious near a phase transition. Coupled 
with the effect of few edges smoothing out the weighting function F(p/h(q)) in the mean field equation, 
we find that theoretically we can see the disappearance of discontinuities in the phase transitions. 

This gives us a crude way to estimate numerically whether a graph will admit a discontinuous or 
continuous phase transition. In Table [I] we perform the procedure above, assuming that F(P) is the 
uniform distribution, and estimate the value of (k) for Erdos-Renyi or scale free graphs for which we 
should see the onset of shocks and discontinuous transitions. A sample of what the adjusted mean field 
curve looks like, near the mean field critical point, with low (k) is also shown in Figure [8| Numerical 
estimates suggest this transition should be fairly sharp and we observed this numerically as well. Because 
we do not have a precise way of determining whether a phase transition has been discontinuous, other 
than simply to observe the size of the change in q at each step forward in p, these results should be taken 
as at most semi-quantitative. However, they do suggest that we have correctly identified the mechanism 
for the disappearance of discontinuous phase transitions. 



graph type 


A = 2 


A = 3 


A = 4 


scale free, 7 = 2.5 


20/20 


14/14 


10/10 


scale free, 7 = 3 


18/20 


12/14 


10/10 


scale free, 7 = 4 


18/20 


12/14 


10/10 


Erdos-Renyi 


18/18 


13/13 


10/10 



Table 1: The theoretical (blue) vs. numerical (black) values for (k) at which we expect to see the 
onset of discontinuous phase transitions (into the q = state) for the uniform distribution. Because these 
predictions require a precise understanding of multistability at the critical point, they should not be taken 
too seriously: we believe a standard deviation of ±2 on each data point is not unreasonable (as we can 
only make scale free graphs with even (k) with our fast algorithm). 
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Figure 8: A sample of how small network effects alter the mean field equation, assuming that F(P) 
corresponds to the uniform distribution, the network is Erdos-Renyi, and h{q) = 1 + 2q. The given value 
of p = 1.16 is quite close to the mean field critical point p c = 1.125. We have added the l/2(k) correction 
to the finite (k) curve. Note that while for the ideal case, the only solution is q = 0, the small (k) curve 
still intersects the line at a positive value of q — i.e., the phase transition has not occurred yet. 

Of course, there remains the question of whether or not these fluctuations actually merge the two 
phases together. This is a question beyond the scope of the simple arguments of this paper. We suggest 
that the answer may be that it depends on the network/disorder/interactions. For example, in the case 
with uniform F(P), there is a very sharp discontinuous phase transition as the value of p is lowered from 
above to below the maximum allowed P v . This transition is hard to remove because the coefficient of the 
l/(k) fluctuations is in fact (as all nodes are in x v = 0), and so network fluctuations will not affect this 
transition. Thus at least for these systems, we conclude that there must be 2 distinct phases, although 
one of the transitions between the two phases may be continuous. 

6. Conclusion 



We have presented the binary decision model and explored its major aspects: understanding from both a 
macroscopic and microscopic view the mean field solutions and fluctuations of basic quantities, and then 
understanding the "glassy" multistability phenomena and the consequences of a spectrum of equilibria. 
The simple form of the decision making, as well as the interactions, allowed us to nearly exactly solve 
this model with very simple, physically motivated arguments. However, there are many questions about 
this model that are still open. Firstly, it is unknown how robust the basic features of this model are to 
modifications: for example, heterogeneity among nodes in the interactions h(q), or "thermal" random 
behavior among the nodes. Secondly, while it is unlikely that the multistability effect is a true spin glass 
effect, it is an interesting question whether or not this model's dynamics at "finite temperature" would 
exhibit aging, another characteristic feature of a glass. Thirdly, an investigation of the model on non- 
random graphs such as hypercubic lattices, or on graphs with many loops, may help to shed some light 
on the nature of multistability, as we mentioned earlier. 

It appears that on heterogeneous networks, the binary decision model has fundamentally different 
behavior from the random field Ising model, despite the similarity in the motivation between the two 
models. One can intuitively see this as follows: for the random field Ising model with coupling J, the 
effective shift to P v is up to k v J, which depends on the degree of node v. Therefore we expect the 
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mean field equations to depend on graph structure, unlike in the binary decision model. Further study 
of qualitative differences in behavior between these two models, as well as modified fiber bundle models 
with the ability to "regenerate" nodes as p is decreased, is a worthwhile direction for further study. 

Finally, it is important to understand the ways in which this model can be tested against empirical 
data. Because we have found that many features of the model are so robust against the detailed network 
structure, and many features of the model only depend primarily on a microscopic quantity a, perhaps 
it is possible to find an elegant way to test this theory against data. Such questions are best saved for a 
future paper. 
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